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ABSTRACT 

Unevenly spaced time series are common in astronomy because of the day-night cy- 
cle, weather conditions, dependence on the source position in the sky, allocated tele- 
scope time, corrupt measurements, for example, or be inherent to the scanning law 
of satellites like Hipparcos and the forthcoming Gaia. This paper aims at improving 
the accuracy of common statistical parameters for the characterization of irregularly 
sampled signals. The uneven representation of time series, often including clumps of 
measurements and gaps with no data, can severely disrupt the values of estimators. 
A weighting scheme adapting to the sampling density and noise level of the signal is 
formulated. Its application to time series from the Hipparcos periodic catalogue led 
to significant improvements in the overall accuracy and precision of the estimators 
with respect to the unweighted counterparts and those weighted by inverse-squared 
uncertainties. Automated classification procedures employing statistical parameters 
weighted by the suggested scheme confirmed the benefits of the improved input at- 
tributes. The classification of eclipsing binaries, Mira, RR Lyrae, Delta Cephei and 
Alpha 2 Canum Venaticorum stars employing exclusively weighted descriptive statis- 
tics achieved an overall accuracy of 92 per cent, about 6 per cent higher than with 
unweighted estimators. 

Key words: methods: data analysis - methods: statistical - stars: variables: general. 



1 INTRODUCTION 

Observed data describe signals entangled with measurement 
properties which depend on the instrumentation, scheduling, 
unexpected events and other factors. Some of the limitations 
resulting from the data acquisition can be addressed statisti- 
cally. For example, in order to characterize the true underly- 
ing signal, the sampling, uncertainties and biases need to be 
accounted for. Pursuing the description of intrinsic signals 
facilitates their interpretation and the comparison of results 
obtained by different experiments. 

A significant number of studies were devoted to the es- 
timation of power spectra and modelling of irregularly sam- 
pled time series (e.g., iCarbonell, Oliver & Ballester||1992 



simple remedy to undesirable effects of irregular sampling 
on the characterization of a signal in terms of statistical pa- 
rameters is presented. Separate works describe corrections 
of biases induced by small sample sizes and Gaussian uncer- 
tainties in the calculation of weighted moments and cumu- 



|Koen||2005| |Vio, Strohmer fc Wamsteker||2000| >. The objec- 
tive of this paper is limited to descriptive statistics, such 
as the mean, variance, higher central moments and robust 
equivalents: they can summarize essential features of sig- 
nals in a few numbers employing straightforward computa- 
tions, which makes them excellent precursors of more de- 
tailed analyses like modelling and classification. Herein, a 
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lants ( Rimoldini|2013a|b[ ). 

Unevenly sampled astronomical time series are common 
in ground- and satellite-based observations, and typically in- 
clude time intervals with clustered and scattered data. For 
exam ple, the sampling laws of surveys such as the Hippar- 
co.n (ESA fl997||Perryman et al.| | l997[ ) and the forthcom- 
ing GaiaF ( |Perryman et al.|2001 1 missions are characterized 
by gaps and clumps of measurements on time-scales much 
greater and smaller than the average sampling interval, re- 
spectivelyrl 

Time series sampled irregularly or with varying can- 
dences may lead to very different estimates of statistical pa- 



1 http://hipparcos.esa.int 

2 http://gaia.esa.int 



3 in the case of the Hipparcos data, sources were typically ob- 
served in sequences of 4 to 6 transits separated by 20 and 108 min 
and repeated every 3 to 5 weeks (Ever et al. 11994 1. 
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rameters for the same signal. For example, if a sinusoid is 
sampled mostly at maximum or minimum, the mean is off- 
set by about the amplitude of the signal, the variance might 
be much smaller than expected (since most points sample 
the same region of the signal), and the distribution appears 
very skewed (the few non-clumped measurements form an 
asymmetric tail in the distribution of measurements). If the 
same clump of data was in proximity of the average level, 
instead, the mean and skewness values could be close to the 
correct value (by serendipity), but the variance would be 
much smaller than the true one. 

Sampling-induced biases do not arise from sparsely sam- 
pled data only and they might manifest independently of the 
number of measurements. While time series with more mea- 
surements are generally associated with better coverage of 
signal features, the importance implicitly assigned to differ- 
ent parts of the signal by unweighted estimators is related 
to the relative frequency of measurements. 

In principle, the most accurate and precise statistical 
parameters could be inferred from the model of a signal. 
Alas, models are often complicated, in the attempt to de- 
scribe the signal features under many circumstances, they 
might require lengthy processing and their accuracy cannot 
be guaranteed in all cases. For example, Fourier series are 
well fitted to model periodic signals, but the description of 
sharp features, like those present in EA-type eclipsing bina- 
ries, requires a high number of harmonics, which can overfit 
smoother parts of the signal and cause unrealistic excursions 
in large intervals with no data ( Dubath et al.|2011 l. 

Linear interpolation is one of the simplest methods to 
approximate a model, assuming the signal features are suffi- 
ciently sampled. This work estimates statistical parameters 
by averaging the linear interpolation of functions of time se- 
ries measurements (depending on the specific estimator). If 
the signal can be recognized in time domain, interpolation 
can be performed in time without requiring further informa- 
tion. If a sparsely sampled signal is primarily mono-periodic, 
it can be interpolated in phase by folding the time series 
with the corresponding period. While the interpolated func- 
tion might include profiles with spiky artefacts, its average 
is rather robust and can be expressed as a weighted mean 
which assigns more relevance to scattered than clumped 
measurements. Statistical parameters weighted by such a 
scheme were tested on data from the Hipparcos periodic cat- 
alogue and led to a significant general improvement in the 
accuracy and precision of estimator values and automated 
classification results, with respect to those obtained with the 
unweighted or error-weighted counterparts. 

This paper is organized as follows. After the definition 
of the notation and terminology in Sec. |2| the description 
of weighted estimators as averages of linear interpolations 
is presented in Sec. [3] with weights defined in time and 
phase, including adaptations for low signal-to-noise regimes 
and small sample sizes. The new weighting scheme is ap- 
plied to time series from the Hipparcos periodic catalogue in 
Sec. p] which includes a comparison of the values of statis- 
tical parameters computed with different weighting schemes 
and their effect on automated classification. The conclusions 
are drawn in Sec. [5] followed by a series of appendices with 
more details on the interpolation of (mono-)periodic signals 
(Appendix Hi), the illustration of modelled light curves to 
verify the accuracy of estimators (Appendix iBp , the defini- 



tions of the statistical parameters employed (Appendix [C| 
and additional scatter plots of estimators (Appendices |D| 
andlEl. 



2 NOTATION AND TERMINOLOGY 

For a set of n measurements x = (xi, X2, •••, x n ), the follow- 
ing quantities are defined: 

(i) The population central moments of order r around the 
mean (i are denoted by (i r and res pective cumulants K2 — (12, 
^3 = (13 and k 4 = (14, — 3(i% (e.g., Stuart & Ord 1969). 

(ii) Sample central moments m r = ^2-_-,Wi(xi — x) r jW 
and respective cumulants k r , where x = X^™=i WiXi/W and 

(iii) Standardized skewness and kurtosis are <?i = ks/k 2 

72 = K4/K2, respectively. 



and g^ = ^4/^2, with population values 71 = ks/k 2 and 



(iv) Noise-unbiased estimates (sometimes called 'denoised' 
for brevity) of central moments and cumulants (Rimoldini 



2013b I are denoted by an asterisk superscript. 



(v) No systematics or instrumental errors are considered 
herein and uncertainties are often referred to as errors, 
(vi) The accuracy of an estimator is related to its distance 
from the true value and thus combines the concepts of bias 
and precision, while the classification accuracy rate is the 
ratio between the number of true positives and the total 
number of sources (of a given type). 

(vii) The precision of an estimator is quantified by its dis- 
persion, while the classification precision rate is the ratio 
between the number of true positives of a given type and 
the total number of sources classified as such a type. 



3 METHOD 

The population values of statistical parameters of a contin- 
uous signal are computed by integrating functions of such a 
signal according to the particular estimator (e.g., leading to 
the areas shaded in blue in Fig. [I]). In case of regular sam- 
pling in fine intervals, the integration is replaced by a sum 
of discrete elements which tends to the population result for 
infinitely small intervals. This limit is not necessarily sat- 
isfied when a continuous signal is sampled irregularly. The 
weighting scheme described in Sec. |3.1| recovers the prop- 
erty that the sum of discrete (weighted) terms approaches 
the result of the continuous function, in the limit of very 
dense (not necessarily uniform) sampling. 

Linear interpolation is a simple method to approxi- 
mate a function with a broken line connecting the data 
points, provided the function is sufficiently sampled in time 
or phase. In this section, such a function is represented by 
the expression of the additive terms of a statistical estima- 
tor. For example, a central moment of order r is defined by 
the average of terms of the form (xt — x) r . Results from 
linear interpolation are equivalent to those from an effec- 
tively infinite regular sampling of the interpolated function. 
In Sec. |3.1| it is shown that the average of a linear interpola- 
tion can be expressed as a weighted mean, which naturally 
removes importance from clumped data (oversampling the 
same region of the signal) in favour of more scattered mea- 
surements (probing extended parts of the signal). Central 
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Figure 1. A sinusoidal signal (blue curve) with mean fi is sam- 
pled unevenly by 10 measurements x; (red circles) and shown in 
the panel on the top-left hand side. The other panels illustrate 
different powers of deviations from the mean of the true signal 
(blue curve) and of the measurements (red circles). In particular, 
(xi — fi) r identifies terms associated with the central moments 
of order r: these estimator values are related to the areas en- 
closed by the blue curves and the zero level (shaded in blue). The 
weights defined by Eqs (|A5j IA7B define the variable widths of 
the shaded bars, which reduce the contribution of clustered mea- 
surements and increase the one of scattered data. In this case, 
although most measurements sample the first half of the signal, 
the increased weights assigned to the remaining data provides a 
better estimate of the area enclosed by the second half of the 
signal than using unweighted schemes. 



moments of a sinusoidal signal are illustrated in Fig. [T] and 
they are related to the areas enclosed by (sin<^) r , which 
are approximated by terms (bar heights) with weights (bar 
widths) adapting to the sampling density in phase. 

Weighting effectively decreases the sample size, since 
more importance is given to some measurements at the ex- 
pense of other ones (e.g., a large sample with only a few 
important elements will be similar to one with only these 
few elements) and exploiting correlations in the data with 
weights (e.g., assigning small weights to measurements sepa- 
rated by small intervals because their values are expected to 
be similar) might introduce biases if applied to expressions 
assuming independent data. However, small biases could 
be justified by significant improvements in precision and a 
mixed weighting scheme which balances precision and accu- 
racy depending on sampling, signal-to-noise (S/N) level and 
sample size is described in Sec. |3.1.3| 



3.1 Weighting schemes 

If the targeted features of signals are sufficiently sampled in 
time, as for the data from the Kepler I Borucki et al.pOlO l 
and CoRoT (Auvergne et al. 20091 missions, linear inter- 



mono-periodic, multi-periodic as well as non-periodic time 
series, and weights can be expressed in time domain. If typ- 
ical time intervals between measurements are larger than 
the resolution needed to sample the signal (as it is often the 
case in the Hipparcos and Gaia surveys) and if the signal 
is mono-periodic or dominated by a single period, sampling 
can be improved by folding the light curve with the value 
of the fundamental period and computing interpolations (or 
weights) in phase. 

The interpolation-based weighting scheme associates 
small and large weights with measurements surrounded by 
small and large data gaps, respectively. If gaps in time or 
phase occasionally cover a relevant fraction of a signal, mea- 
surements at the gap boundaries cannot provide sufficient 
information on the features of the signal and the impor- 
tance of such data could be limited, for example, to some 
maximum weight value. This strategy can extend the appli- 
cability of the following weighting scheme to under-sampled 
signals, periodic or non-periodic, to mitigate statistical bi- 
ases due to data clumps on scales smaller than the signal 
features. 



3.1.1 Weights in time 

If an estimator is defined by the average of a function 
in time, it can be approximated by the mean of the linear 
interpolation of terms Oi at times ti, where i G (l,n) for a 
time series of n measurements sorted in time, and it can be 
expressed as a weighted average as follows: 



1 ^ /"**+! 



ti+1 — t, 



(t - ti) 



dt 



+ Oi+l 



Zn t 



_ £ - '-^ {U + 1 _ U ) 



2(t n -tl) 



2(t„-tl) 



/ J 0i (ti+l — ti) + 2_^ 0i (ti — ti-1 

i—1 i—2 

n-1 

y J 0i (ti+i — ti-i) + 

1=2 

+ 01 (t'2 — tl) + n (tn — tn-l) 



(1) 

(2) 
(3) 



1 " 



where 

Wi — t i+ i — ti-i Vi € (2, n — 1) 
Wi — t-i, — t\ 

Wn — tn tn — 1, 

and 

n 

W / = ^«) l = 2(i n -i 1 ). 



(4) 
(5) 



(6) 
(7) 
(8) 



(9) 



polation is expected to provide reliable approximations of 



When differences between successive times ti are too large 
for sensible interpolation, they might be limited to some 
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maximum interval At n 



as follows: 



Wi — min {U+i—U, Ai max } + 

+ min{ti-ti-i,At ma =} Vie(2,n-1) (10) 

wi —min{t2 — h, At max } (11) 



W n — nim 



{*„ 



n-l 



^t m axj 



(12) 



and, of course, then W ^ 2(t n — £i). The maximum interval 
Ai max might be expressed as a function of n to better control 
the precision of the estimator in the limits of both high and 
low values of n, as follows: 



At max = h(n\a,b) 

where 

h(n\a, b) 



for a, b > 0. 



(13) 



(14) 



1 + e~( n - a )/ b 

The function h(n\a, b) constitutes just an example to achieve 
a mixed weighting scheme: tuning parameters a, b regulate 
the transition from unweighted to weighted estimators in the 
limits of small and large n, respectively. The determination 
of such parameters might involve the maximization of the 
overall accuracy of estimators employing data simulated in 
a context similar to the one intended for analysis. 

3.1.2 Weights in phase 

If sampling in time is sparse but the data exhibit periodicity, 
the time series can be folded with the dominant period to in- 
crease the average sampling rate in phase (by a factor equal 
to the duration of the time series divided by the period) . 

Since interpolation of phase-sorted data is carried over 
from the last to the first point in phase, weights in phase <j> 
corresponding to Eqs ( |10[ )-( [T2| become (see Appendix |A|| : 

Wi — min {<j>i+i — fa, A0 max } + 

+ min{(/>i — <jh-i, A0 max } Vie(2,n— 1) (15) 

lUi = min {(f>2 — 4>i, A0 max } + 

+ min {0i - <f) n + 2n, A0 max } (16) 

W n =min{0„ — 4> n -i, A0 max } + 

+ min {0i — <j> n + 27r, A0 max } (17) 

with W ^ 4-7T and, for example, A0 max = 2ir h(n\a' , b'). 

3.1.3 Weights with noise and sample-size dependence 

In order to avoid interpolating large noise fluctuations, 
weights could be set to reduce to inverse-squared uncertain- 
ties at low S/N, as they proved more precise in the simu- 
lations described by Rimoldini (2013b I. The transition be- 



tween high and low S/N regimes could be pursued with a 
function h(S/N\a,b) as defined in Eq. (14 1, depending on 



the S/N ratio and tuning parameters a, b. Another function 
of the form h(n\a' , b') could be used to reduce the relevance 
of interpolation-based weights for small sample sizes n, as 
follows: 



Wi = h(S/N\a,b) 



« ! < 



Era 
.7=1' 



+ [l-h{S/N\a,b)] 



£• 



Wj — h(n\a ,b ) 



e; 



-+[l-h(n\a',b')]/n, 



(18) 
(19) 



where e* denotes the uncertainty of the i-th measurement 
and w'( is defined by Eqs (|6|-(l8l or (A5l-(A7l. Tuning pa- 
rameters a,a',b,b' offer the possibility to reach a compro- 
mise solution between precision and accuracy at high and 
low values of n and S/N ratios, according to the specific 
estimators, signals, sampling, errors, sample sizes and their 
distributions in the data. An application of such tuning pa- 
rameters to simulated data, as a function of sample size and 
S/N ratio, was illustrated in Rimoldini (2013a|b'l. 



4 HIPPARCOS PERIODIC VARIABLE STARS 

The effect of interpolation-based weights on statistical pa- 
rameters was explored for a realistic distribution of signals, 
represented by time series from the Hipparcos catalogue of 
periodic variable stars. The Hipparcos mission (ESA 1997) 
performed astrometric and photometric measurements of 
the brightest sources in the sky. The full catalogue (|Perry-| 



man et al.|1997l contains 118 204 sources with photometry. 



Among the 11 597 stars identified as variables, a reliable pe- 
riod could be computed (or was consistent with the one from 
literature) for 2 712 objects, which were published in the pe- 
riodic catalogue (Vol. 11 |ESA|1997| |Eyer|1998| ). This set of 
sources, which contained light variations dominated mostly 
by single periods, was chosen to illustrate the application of 
one of the weighting schemes described in Sec. [3] 



4.1 Statistical parameters 

In order to assess the accuracy of estimators with respect 
to the (unknown) real signal, the latter was assumed to be 
represented by a model of the time series and data were 
simulated employing measurement uncertainties as Gaus- 
sian random variables (around the model) and phases as 
in the real data. Only good quality measurements were ac- 
cepted, flagged by the field HT4 as zero or one (see Vol. 1 of 
|ESA|1997[ ) , which reduced the number of sources considered 
to 26830 



4-1.1 Light-curve models 

Time series were folded with the period provided by the Hip- 
parcos cataloguajand modelled by a cubic smoothing spline 
(the smooth, spline function from the stats package in R, 
R Development Core Team|2013 |. The condition of period- 
icity at the boundaries was approximated by replicating a 
whole cycle of folded data before and after the cycle con- 
sidered as reference for the 'true' statistical parameters and 
for the generation of simulated data. The smoothing param- 
eter was estimated from the data with a generalized cross- 
validation method (e.g., Ruppert, Wand, & Carroll 20031 
and the degrees of freedom were adjusted in special cases, 



4 The Hipparcos identifiers of sources with no good quality mea- 
surements (i.e., with field HT4 > 1 only) were: 1196, 10027, 17878, 
20570, 24019, 25673, 39084, 42715, 42726, 46502, 52538, 53937, 
58112, 60904, 61997, 63125, 69582, 72583, 88905, 90026, 93595, 
93724, 96007, 99675, 102246, 10240 9, 112317, 1 12470 and 118188. 

5 The mean period (field Pll of [ESA][l997) derived from the 
Hipparcos data was employed. If this was not available, the period 
from the literature (listed in field P18) was considered. 
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to mitigate over-fitting highly clumped data with large gaps 
as well as under-fitting the profile of eclipsing binariesfjThe 
combination of large gaps and sharp features did not lead to 
accurate models and other methods might help avoid mod- 
elling artefacts in gapped data, although these are out of 
the scope of this application. Also, cases in which sparse 
sampling missed important features could not be improved. 
Nevertheless, the smoothed best-fitting curves seemed to 
capture the relevant shapes of true signals in most cases. 

The resulting models are presented together with the 
original and simulated data in Appendix [B] The difficulty 
to achieve accurate models for all sources and the need of 
more complex modelling techniques were confirmed. While 
models were sometimes not ideal, most of them were of suffi- 
cient quality to supply a realistic distribution of the relevant 
features of the Hipparcos periodic variable stars. Differences 
from data of other surveys were expected to be greater than 
those due to modelling inaccuracies. Less than one per cent 
of all sources (23 time seriesfj with significant modelling 
artefacts in gaps with no data were removed (many of them 
were EA-type eclipsing binaries with large data gaps), thus 
2660 objects were included in the assessment of statistical 
parameters. The number of measurements per time series 
ranged from 18 to 331, with S/N ratios from 0.2 to 116, 
according to the definition in Eq. (1201). 



4-1-2 Weighting scheme 

The Hipparcos light curves were folded with the catalogue 
period and interpolation-based weights were computed in 
phase. For brevity, such weights are referred to as 'phase 
weights', while inverse-squared uncertainties are called 'error 
weights'. 

The estimators employed herein included weighted mo- 
ments and cumulants corrected for biases from Gaussian un- 
certainties (called 'noise- unbiased' or 'denoised' estimators; 
see |Rimoldini|2013b[ ) and some robust weighted measures, as 
defined in Appendix [C] Simulations in Rimoldini ( 2013a|b l 
suggested that noise-unbiased phase-weighted sample mo- 
ments can be more accurate for S/N > 2 and sample sizes 
7i > 20 with respect to other schemes, while error weight- 
ing appeared the most appropriate option for noisy signals. 
Thus, the weighting scheme chosen for this application com- 
bined error and phase information as described by Eqs ( 18 1 



6 The number of degrees of freedom df depended on the un- 
weighted standardized sample skewness S of the data, which iden- 
tified light variations typical of eclipsing binaries, and the ratio 
Ti between the median and the third largest gaps in phase, to 
better deal with clumped data with large gaps. Denoting by n 
the number of measurements in a time series, df equalled 15, 24 
and 36 for 1Z less than 6/n, 10/n and 20/n, respectively, provided 
S < 1.5. For greater values of 5, df was set to the smallest integer 
not less than 3ra/5. The quoted values of df take into account the 
replicated data at each of the extremes of the folded light curve 
(to induce quasi-periodic boundary conditions on the model of 
the central cycle). 

7 The Hipparcos identifiers of the 23 sources removed from con- 
sideration because of poor modelling were: 1901, 4279, 21600, 
23416, 23453, 25591, 32397, 40853, 42853, 45094, 48054, 58854, 
61281, 68064, 73533, 76152, 89579, 90313, 95611, 96739, 104483, 
108317 and 112928. 



and ( A5 i-( A7l. The balance between phase weights (at high 
S/N) and error weights (at low S/N) was controlled by the 
parameters a = 3 and b = 1.2, which provided a satisfactory 
overall accuracy for the set of estimators and light-curve 
shapes of the Hipparcos periodic variable stars. The signal- 
to-noise ratio S/N was estimated from the true (model) sig- 
nal variance /Z2 and the average of squared measurement 
uncertainties ei as follows: 



S/N = 



l'2 



-,1/2 



^Efci- 



(20) 



4-1.3 Results 



The deviations of error-weighted sample statistical parame- 
ters from population values (computed from a fine regular 
sampling of the models) are compared to the noise-unbiased 
error-phase weighted counterparts as a function of S/N ratio 
in Figs[2|j3] Error-phase weighted (and noise-unbiased, when 
applicable) estimators were generally more accurate than 
error-weighted estimators. Since larger light variations were 
correlated with higher S/N levels, the effect of the correla- 



tion between errors and magnitudes (e.g., see van Leeuwen 



19971 was visible when weighting by errors at high S/N: the 
mean was biased towards brighter measurements, the vari- 
ance tended to be smaller than the real value and the scatter 
of higher moments around the true values was greater. Such 
effects were significantly alleviated by error-phase weights. 
Also, the accuracy of the noise-unbiased variance m-j and 
kurtosis moment m| were much improved at S/N < 3. In 
the particular case of the standardized kurtosis, the value of 
mi/mi — 3 — 72 equalled the one of the cumulant k±/m\ — 72 
by definition, for any weighting scheme (unlike the noise- 
unbiased counterparts) r] 

Estimators which provide similar information (such as 
the mean and median, variance and interquartile range, 
skewness and kurtosis standardized by the estimated or true 
variance) are compared in scatter plots in Appendix [D] Al- 
though contributions from different S/N ratios are not dis- 
tinguished, error-phase weighted (and noise-unbiased, when 
applicable) estimators are always associated with more 
strongly peaked distributions around the true values than 
the error-weighted analogues. 

In order to quantify the effect of error-phase weighting 
and additional denoising with respect to simple error weight- 
ing, the top panels in Figs[4j|6]indicate the fraction of sources 
associated with estimators improved by phase weights and 
noise correction. The accuracy of a generic estimator £ is as- 
sessed by its distance \A£ | = \£ — 77 1 from the true value n. 
The abscissas in Figs HJjSl represent the difference between 
the distance of error-phase weighted (and noise-unbiased, 
when applicable) estimators and the one of error-weighted 
estimators from the correct values: negative differences indi- 
cate smaller distances to the true values and thus improved 
accuracy with respect to error-weighted estimators. The ac- 
curacy of error-phase weighted and noise-unbiased estima- 
tors improved in 70-to-90 per cent of the cases and deteri- 



The deviations from model values of the standardized noise- 



ization by (m?;) 2 . 
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Figure 2. Deviations of the mean, median, variance, interquartile range (IQR) and skewness from their true values, as a function of S/N 
ratio. Green triangles denote error-weighted estimators, red crosses indicate error-phase weighted estimators and blue circles represent 
noise-unbiased error-phase weighted estimators. In the case of the median and IQR, 'true' refers to the true median and IQR values, 
respectively. 



orations were typically smaller in frequency and magnitude 
than improvements. In the singular case of the standardized 
skewness g\ , denoising worsened almost 8 per cent of the es- 
timates with respect to error-phase weighting. As apparent 
in Fig. [5k, the accuracy of the standardized noise-unbiased 
skewness worsened at low S/N ratios. This was expected 
because the noise-unbiased variance might be easily under- 
estimated when uncertainties are of the same order as the 
signal, and the normalization of the skewness by a much 
smaller number than the correct one could rapidly degrade 
the precision of the estimator. A similar (less pronounced) 
tendency at low S/N ratios was also observed in the case of 
the standardized kurtosis. 



As noted in Rimoldini (2013b I, the noise-unbiased vari- 



ance can become too small at low S/N levels, overestimating 



the standardized noise-unbiased skewness and kurtosis, or 
leading to an undetermined skewness value in case of non- 
positive variance. On the other hand, when the variance is 
not corrected by noise biases, it is generally overestimated, 
and if the true values of skewness or kurtosis are sufficiently 
close to zero, the standardized noise- biased skewness and 
kurtosis become more accurate (by serendipity). 



4.2 Automated classification 

The effect of statistical parameters weighted by differ- 
ent schemes on automated classification was assessed by 
comparing the classification accuracy and precision as a 
function of variability type employing unweighted, error- 
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Figure 3. Deviations of the kurtosis moment and cumulant from their true values, as a function of S/N ratio. Green triangles denote 
error-weighted estimators, red crosses indicate error-phase weighted estimators and blue circles represent noise- unbiased error-phase 
weighted estimators. 



weighted, phase-error weighted and noise-unbiased phase- 
error weighted estimators. 



4-2.1 Attributes and variability types 

Automated classification of stellar variability types was pur- 
sued with a set of attributes which characterized features of 
different classes. Some studies employed only information 



from light-curve modelling (Debosscher et al. 2007 2009 



Blomme et al. 2010 20111, with additional statistical pa- 



rameters ([Richards et al.|2011 1 or colour information ( Sarro 



et al.|2009 1 , while attributes from modelling, statistical and 



astrophysical quantities were used in Dubath et al. (20111; 



Rimoldini et al. ( 2012 1 



The list of attributes employed for classification herein 
was restricted to descriptive statistics (mean, median, vari- 
ance, interquartile range, normalized and non-normalized 
skewness and kurtosis), in order to test the effect of dif- 
ferent weighting schemes on classification. Only variability 
types which could be identified by the distribution of mea- 
surements were included in the training set, such as eclipsing 
binaries, RR Lyrae, Mira, S Cephei and Alpha 2 Canum Ve- 
naticorum stars, accounting for over 60 per cent of all sources 
in the Hipparcos periodic catalogue. 

In these classification experiments, statistical parame- 
ters did not depend on modelling and were computed on 
the original (not simulated) data, which thus included the 
sources previously excluded because of modelling issues. The 
same quality flags, periods and weighting schemes were ap- 
plied as described in Sec.|4.1| Since the true signal variance 



was unknown, the S/N ratio was estimated by substituting 
the unknown /12 in Eq. ( |20[ ) with the noise-unbiased phase- 
weighted sample variance m| and by weighting the average 
of squared errors in the denominator with the same weights 
for consistency: 



S/N 



J27=l W i( X i 



+ Efei«?e?/W r 



En 9 



1 



1/2 



(21) 



The S/N ratios of the sources selected for classification (as 
described in Sec. 4.2.21 spanned a range from 0.4 to 140. 



4-2.2 Data selection 

Sources from the Hipparcos periodic catalogue were cross- 
matched with classificati ons from literature as available i n 
the Variable Star Index ( |Watson, Henden fc Price||2012p 
of the American Association of Variable Star Observers 
(AAVSO) with the nearest object within 1 arcsecrj 

Uncertain classifications (with class labels followed by 
the mark ':') were excluded, unless the uncertainty referred 



9 See http://www.aavso.org/vsx/index.php?view=about.top for 
more details on how literature information has been selected, 
maintained and revised. A comprehensive list of variability 
types, labels and their descriptions is available at http://www. 
aavso.org/vsx/help/VariableStarTypeDesignationsInVSX.pdf 
"' Only one source, Hip 31400, was not considered because asso- 
ciated with different classes at very similar angles from the direc- 
tion of the Hipparcos source. 
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Figure 4. The absolute deviations from population values of error-phase weighted estimators (red crosses and dashed lines, or black 
circles and black solid lines) with additional denoising, when applicable (blue circles and blue solid lines), are denoted by a prime symbol 
and their difference from the error-weighted counterparts is presented as a function of S/N ratio. Negative values indicate improvements 
in accuracy with respect to error- weighted estimators. The fraction and magnitude of improved cases can be inferred from the cumulative 
distributions in the top panels. The following estimators are included: (a) mean, (b) median, (c) variance and (d) interquartile range. 



to properties of eclipsing binaries other than light-curve 
shapes, such as the physical characteristics, the luminos- 
ity class of the components, or the degree of filling of the 
inner Roche lobes. Objects associated with combinations of 
different classes (with labels joined by the symbol '+') were 
not addressed herein. The subset of variability types chosen 
to test classification with statistical parameters are listed in 
Tableland include 1605 sources associated with class labels 
EA, EB, EW, ACV, M, RRAB, RRC, DCEP and DCEPS, 
defined in Table [l] together with the corresponding sample 



sizes. Sources which poorly represented their class were not 
removed from the training set in order to avoid the intro- 
duction of selection biases (e.g., by favouring objects with 
higher S/N ratios) and because the focus of this application 
was related to the relative classification accuracy employing 
attributes with different weighting schemes. 
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4-2.3 Training-set attributes 

The distributions of a selection of error-phase weighted (and 
noise-unbiased, when applicable) statistical parameters as 
a function of variability type are presented in Fig. [7] The 
mean or median magnitudes of sources at broadly different 
distances from the observer do not provide information on 
the intrinsic source properties. However, the correlation be- 
tween magnitude and noise, coupled with the interquartile 
range or variance (as shown in Fig.lTa), is related to the S/N 
level and thus to an observational selection. Figure [7b illus- 



trates clearly the difference between robust and non-robust 
estimators for the same quantity, such as the interquartile 
range versus the variance, which proves effective at separat- 
ing eclipsing binaries (mostly EA and some fraction of EB 
types). The interquartile range and standardized skewness 
are presented in Fig.[7fc: the skewness separates the eclipsing 
binary subtypes from most other classes, which are better 
distinguished by the interquartile range. A similar scatter 
plot is shown in Fig. [TH in terms of normalized and non- 
normalized kurtosis moments, with the difference that some 
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Table 1. The training set employed to test the effect of descrip- 
tive statistics weighted by different schemes on classification in- 
cludes 9 variability types and 1605 sources from the Hipparcos 
periodic catalogue. 



Variability Type 



Label 



Number 



Eclipsing Binary: Algol type 

Beta Lyrae type 

W Ursae Majoris type 
Alpha 2 Canum Venaticorum 
Mira Ceti 
RR Lyrae: Asymmetric light curve 

Nearly symmetric light curve 
Delta Cephei 

First overtone pulsators 



409 


208 


170 


183 


218 


139 


28 


220 


30 



Total: 



1605 



classes occupy different relative loci (such as the EA types 
and the RRAB with respect to DCEP variables). 

The estimators and variability types illustrated in Fig. [jj 
are also presented for the unweighted case in Appendix [Ej 
The distributions of most unweighted estimators in Fig. |E1| 
tend to be more scattered than the ones shown in Fig. [jj In 
the case of ACV stars, instead, the unweighted skewness and 
kurtosis are less scattered, because these sources have very 
low S/N ratios (typically S/N < 2.5) and the correction 
of noise biases at such S/N levels can decrease significantly 



the precision of higher moments ( Rimoldini|2013b I . Another 
noticeable difference between Figs |7| and |E1| is related to 
the larger skewness g\ of Mira stars in the unweighted case, 
which can be understood by the systematic difference in the 



uncertainties associated with faint and bright measurements 
(e.g., [van Leeuw cn 1997). Such a difference is enhanced by 
the large amplitudes of light variations of Mira types, so that 
greater uncertainties on the faint side of the signal gener- 
ate more faint than bright 'outliers', biasing the unweighted 
skewness towards greater values. 



4-2.4 Classification method 

Automated clas sification tests were performed employing 
random forests ( Breiman||200l| ) | n | Random forest is an ac- 
curate tree-based classification method (e.g., see |Richards| 
et al.|2011 l, quite robust to outliers and strongly correlated 
attributes. The accuracy of the classifier was estimated from 
a subset of sources (about one-third of all objects) randomly 
omitted from the learning process and thus called 'out-of- 
bag' sources. The importance of an attribute for classifica- 
tion was measured by the mean decrease in accuracy af- 
ter permuting the values of that attribute in the out-of-bag 
sources (see 



Rimoldini et al.|2012 for more details) 



The random forest classifier was trained with estima- 
tors employing a single weighting scheme per run, i.e., 
unweighted, error-weighted, error-phase weighted or noise- 
unbiased error-phase weighted. For each weighting scheme, 
random forest was executed 1000 times with 500 trees and 
the classification accuracy and precision rates were aggre- 
gated after each run in order to assess their mean and dis- 
persion. 



This work employed the randomForest package implemented 



in R (R Development Core Team 2013 I by 



Liaw & Wiener 



(2002) 
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Figure 7. A selection of estimators employed for classification is illustrated as a function of variability type. All of the estimators are 
phase-error weighted (primed) and some are also noise-unbiased (starred), as defined in AppendixO Class labels are described in Tablcll] 
and denoted by symbols as shown in the legend of each panel. 



4-2.5 Results 

The mean classification accuracy and precision rates are 
listed as a function of variability type and attribute weight- 
ing scheme in Tables [2] and [3] and illustrated in Figs [8k 
and [8b, respectively. Training-set sources were reclassified 
by sole statistical parameters with an overall accuracy from 
77.0 to 82.6 per cent on average, depending on the weighting 
scheme adopted. Unweighted estimators led to the smallest 
accuracy level, closely followed by error-weighted estimators 
at 77.6 per cent. The introduction of phase weights improved 
the accuracy to 82.3 per cent and noise-unbiased error-phase 
weighted estimators achieved the best overall classification 



accuracy of 82.6 per cent (with uncertainty at the level of 0.2 
to 0.3 per cent). The overall trends were generally reflected 
by single variability types with only a few exceptions. In 
particular, error-weighted estimators of RRAB types led to 
significantly worse classification accuracy than unweighted 
estimators. This was explained by the relatively large am- 
plitude of RRAB variables, which enhanced the systematic 
difference of uncertainties between bright and faint measure- 
ments in the light curve: weighting by errors decreased the 
importance of faint measurements, increasing the similarity 
(and thus confusion) with DCEP and M-type light curves. In 
the case of DCEPS variables, denoising was marginally coun- 
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Table 2. Classification accuracy rates (per cent values) from the 
reclassification of training-set sources from the Hipparcos peri- 
odic catalogue are listed as a function of variability type and 
attribute weighting scheme. Class labels are defined in Table ffl 
Accuracies and corresponding uncertainties represent average val- 
ues from 1000 runs of random forests employing 500 trees. 



Var. 




Error 


Error- Phase 


Err.-Ph. w. 


Type 


Unweighted 


Weighted 


Weighted 


& Denoised 


EA 


87.7 ± 0.3 


87.7 ± 0.4 


88.9 ± 0.4 


89.0 ± 0.3 


EB 


53.7 ± 1.1 


56.4 ± 1.0 


58.7 ± 0.9 


59.8 ± 1.0 


EW 


59.2 ± 1.4 


60.9 ± 1.2 


68.5 ± 1.0 


70.2 ± 1.0 


ACV 


84.9 ± 0.9 


85.9 ± 0.9 


88.8 ± 1.0 


87.6 ± 0.7 


M 


96.7 ± 0.2 


96.5 ± 0.3 


98.0 ± 0.3 


98.2 ± 0.3 


RRAB 


83.1 ± 0.8 


79.4 ± 1.0 


88.3 ± 0.6 


87.8 ± 0.8 


RRC 


41.2 ± 3.0 


47.0 ± 4.0 


41.8 ± 3.1 


50.8 ± 3.3 


DCEP 


73.9 ± 0.8 


74.7 ± 0.7 


86.6 ± 0.6 


87.1 ± 0.6 


DCEPS 


31.3 ± 3.5 


35.1 ± 3.8 


56.6 ± 2.2 


51.4 ± 1.9 


ALL 


77.0 ± 0.3 


77.6 ± 0.3 


82.3 ± 0.2 


82.6 ± 0.2 



Table 3. Classification precision rates (per cent values) from the 
reclassification of training-set sources from the Hipparcos peri- 
odic catalogue are listed as a function of variability type and 
attribute weighting scheme. Class labels are defined in Table [l] 
Precisions and corresponding uncertainties represent average val- 
ues from 1000 runs of random forests employing 500 trees. 



Var. 




Error 


Error- Phase 


Err.-Ph. w. 


Type 


Unweighted 


Weighted 


Weighted 


& Denoised 


EA 


88.5 ± 0.4 


89.4 ± 0.4 


91.2 ± 0.4 


91.0 ± 0.3 


EB 


58.4 ± 1.0 


59.9 ± 1.0 


63.1 ± 0.9 


64.5 ± 0.8 


EW 


58.5 ± 1.1 


59.4 ± 1.0 


65.5 ± 0.9 


65.2 ± 0.9 


ACV 


72.4 ± 0.5 


74.8 ± 0.6 


76.3 ± 0.5 


76.6 ± 0.5 


M 


98.2 ± 0.2 


97.2 ± 0.2 


98.6 ± 0.3 


98.8 ± 0.3 


RRAB 


80.8 ± 0.7 


80.5 ± 0.9 


89.8 ± 0.8 


90.5 ± 0.9 


RRC 


52.5 ± 3.2 


56.5 ± 3.3 


82.4 ± 4.9 


88.4 ± 4.3 


DCEP 


74.7 ± 0.8 


74.7 ± 0.8 


85.3 ± 0.6 


85.4 ± 0.6 


DCEPS 


37.9 ± 3.5 


39.7 ± 3.3 


55.4 ± 1.9 


56.2 ± 1.8 


ALL 


77.0 ± 0.3 


77.6 ± 0.3 


82.3 ± 0.2 


82.6 ± 0.2 



terproductive for accuracy, although still significantly better 
with respect to the unweighted or error-weighted schemes. 
This was expected to be related to the increased confusion 
with the noise-unbiased estimators of EW types, which were 
more numerous and extended to lower S/N ratios than the 
DCEPS stars: while the distributions of the standardized 
error-phase weighted kurtosis moments of DCEPS and EW 
types could just be separated, the additional denoising de- 
creased the values of standardized kurtosis moments more 
for the EW than the DCEPS stars, leading to overlapping 
distributions dominated by the EW types. 

Precision rates of classification per variability type 
showed successive improvements (consistent within uncer- 
tainties) from the unweighted to the error-weighted, error- 
phase weighted and noise-unbiased schemes. A single excep- 
tion was related to Mira stars: their classification precision 
with error-weighted estimators was about one per cent worse 
than with unweighted attributes, as a consequence of the in- 
creased contamination by RRAB stars (as explained above). 
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Figure 8. Classification accuracy and precision rates are shown 
in panels (a) and (b), respectively, as a function of variability type 
and attribute weighting scheme. The values are obtained from the 
reclassification of training-set sources from the Hipparcos periodic 
catalogue and are listed in Tables[2]and|3] Class labels are defined 
in Tablellland colour coded as shown in the legend of each panel. 



The importance of phase- weighted (and noise- unbiased, 
when applicable) attributes as a function of variability type 
from a run of random forest is depicted in Fig. [9] with darker 
cells indicating more important attributes for a given class. 
For example, the skewness confirmed to be particularly use- 
ful to distinguish eclipsing binaries and the asymmetric light 
curves of RRAB stars. Also, large and small amplitude vari- 
ables, such as ACV and M types, could easily be separated 
by the variance. 

The confusion matrix from a run of random forest, em- 
ploying error-phase weighted (and noise-unbiased, when ap- 
plicable) estimators, is shown in Fig. |10| The overall classi- 
fication accuracy was 83 per cent and reached 92 per cent 
after aggregating the subtypes of eclipsing binary, RR Lyrae 
and Delta Cephei into their superclasses. The improvement 
in classification accuracy, with respect to results employing 
unweighted estimators, was about 6 per cent in both cases 
of separated and aggregated variability subtypes. The level 
of confusion between eclipsing binaries of EA, EB and EW 
types was expected from the natural overlap in their defini- 
tions. The misclassification of many ACV types as eclipsing 
binaries was attributed to the restricted set of attributes 
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Figure 9. The importance of attributes is illustrated as a func- 
tion of variability type, with darker grey levels indicating more 
important attributes for the identification of a given class. At- 
tributes are phase-error weighted (primed) and some also noise- 
unbiased (starred), as defined in Appendix [C] Class labels are 
described in Table [T] 



employed herein (e.g., the attribute 'QSO variations' in Ri- 
moldini et al. 2012 proved to be useful in this context) 



Similarly, the confusion between RR Lyrae and Delta Cephei 
variables would not have occurred if the period was included 
in the set of classification attributes. On the other hand, sim- 
ple amplitude estimators such as variance and interquartile 
range were sufficient to separate Mira stars from most other 
classes with high accuracy (as apparent in Fig. Us. 



5 CONCLUSIONS 

A simple weighting scheme based on linear interpolation, 
applicable in phase or time, was shown to improve the accu- 
racy of descriptive statistics of unevenly sampled time series. 
Weights represented by time intervals could be applied to 
well-sampled signals (mono-periodic, multi-periodic, semi- 
periodic and non-periodic). For sparsely sampled but peri- 
odic time series (dominated by a single period), light curves 
could be folded with the fundamental periods and weights 
expressed in terms of phase intervals. 

The Hipparcos catalogue of periodic variables repre- 
sented a suitable test bed of unevenly sampled light curves 
with realistic distributions of variability types to investigate 
the accuracy of estimators weighted by different schemes. 
Noise-unbiased estimators weighted by phase intervals (or 
inverse-squared uncertainties at low S/N ratios) improved 
the accuracy for 70-to-90 per cent of the sources with respect 
to the values of the error-weighted counterparts. Deteriora- 
tions in accuracy were observed in 10-to-30 per cent of the 
cases, with magnitudes typically smaller than the ones cor- 
responding to improvements. 

Automated classification experiments confirmed that 



Figure 10. The confusion matrix of training-set sources reclassi- 
fied by random forests with 500 trees and employing error-phase 
weighted (and noise-unbiased, when applicable) estimators. The 
overall classification accuracy is 83 per cent and reaches 92 per 
cent after aggregating the subtypes of eclipsing binary, RR Lyrae 
and Delta Cephei into their superclasses. Rows and columns refer 
to literature and predicted types, respectively. Class labels are 
defined in Table [J 



the best overall results were achieved employing the set of 
the most accurate attributes, i.e., the noise-unbiased esti- 
mators weighted by phase intervals and inverse-squared un- 
certainties at high and low S/N ratios, respectively. The 
overall improvement in classification accuracy with respect 
to the result employing unweighted estimators was about 6 
per cent and most of it was related to the introduction of 
phase weights which adapted to the varying sampling den- 
sity of the light curve. 
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APPENDIX A: PERIODIC INTERPOLATION 

Substituting time-sorted data with phase-sorted measure- 
ments in Eq. <TlJ> and adding the interpolation term from 
the last to the first point in phase, it follows: 



2tt 

1 

4-k 

1 

4-k 



1 v^ 0i + 9 i+1 

-^ ;l — 2 — (0i+i _ *° 

-i 

n n+1 

y^ Oi (<Pi+i - <Pi) + 2^ Oi {fa - 4>i- 

_i=l i=2 

n 

> y ft (4>i + l — 4>i-l) + 
_i=2 

+ 6l (02 — 0l) + fti+1 (0n + l — 4>r, 

1 - 

i — l 

where # n +i = #i and n +i — 01 + 27r, so that 
Wi — (j> i+ i — <f)i-\ Vi G (2, n — 1) 

w n — 01 — 0„_1 + 27T, 

and 

71 

W — y~2 Wi = 4vr - 



(Al) 
(A2) 



(A3) 
(A4) 



(A5) 
(A6) 
(A7) 



(A8) 



When differences between successive phases are a significant 
fraction of a cycle, they might be limited to some maximum 
interval A0 max as follows: 

Wi — min {<f>i+i — (pi, A0 max } + 

+ mm{<pi-^i- u A0 max } Vi€(2,n-1) (A9) 
wi — min {02 — 0i, A0 max } + 

+ min {0! - 0„ + 2tt, A0 max } (A10) 

W n =min{0 n - 0„_i, A0 max } + 

+ min {0i -(/>„ + 2tt, A0 max } (All) 

and, of course, then W ^ 4n. 



APPENDIX B: THE HIPPARCOS PERIODIC 
LIGHT CURVES 

Light curves of the Hipparcos periodic variables are illus- 
trated in Fig. |B1| Panels on the left-hand side show the 
Hipparcos data folded with the catalogue period, while the 
right-hand panels present data simulated around the model 
according to the measurement uncertainties (assumed Gaus- 
sian) and the observed phases. Simulated data were em- 
ployed to assess the accuracy of statistical estimators (with 
respect to the model), while classification was performed on 
the original data. Figure [rJl] provides sample light curves for 
4 of the 2683 Hipparcos sources available online (see sup- 
porting information)! 12 1 



12 The 29 objects with no data associated with good quality flags 
(i.e., field HT4 > 1 only) were not included in Fig 
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Figure Bl. Panels on the left-hand side present folded light curves from the Hipparcos periodic catalogue, together with the smoothing 
spline model employed to generate simulated data (shown in the panels on the right-hand side). Simulations served to assess the accuracy 
of statistical parameters with respect to the model (as reference), while classification experiments were performed on the original data. 
The illustration of the full set of light curves employed herein is available online (see supporting information). 
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APPENDIX C: STATISTICAL PARAMETERS 

The statistical estimators employed to test different weight- 
ing schemes are defined below. 



CI Variance, skewness and kurtosis 

Sample noise-unbiased weighted moments and cumulants 



are defined following Rimoldini (2013b I. Denoting the sam- 



ple weighted central moments of order r by 
1 - 

i=l 

where 

n 
i = l 

and representing the Gaussian uncertainties by et, the noise- 
unbiased estimators (starred) are defined as follows: 
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(C6) 

(C7) 
(C8) 
(C9) 



fc 4 =m\ — 3 (m 2 )* 

ff i=fcl/(fc 2 *) 3/2 

C2 Percentiles 

Percentiles depend on the rank of sorted values, thus they 
are less sensitive to extreme values than moments and cu- 
mulants which involve powers of deviations from the mean 
and average over all elements. The m-th percentile P m (x) 
is defined as the (interpolated) value such that m per 
cent of the data are smaller than P m (x). Two common 
percentiles are the median Pso(x) and interquartile range 
IQR = P 75 (x)-P2 5 (x). 

Denoting the list of measurements Xi sorted in in- 
creasing values by {k(i), .. .,#(„)}, associated with weights 
{W(i), ..., U>(n.)}) respectively, the m-th weighted percentile 
P m (x) is defined as follows: 



Pm(x) = 




C(i+1) 



Hi)) 



if < m < pi 
if pi < m < p,+i 
if p„ < m < 100 
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Figure D2. Deviations from the true values of the kurtosis cumu- 
lants, standardized by the estimated and true variances. The tri- 
angles and histograms in green denote error- weighted estimators, 
the crosses and histograms in red indicate error-phase weighted 
estimators, and the circles and histograms in blue represent noise- 
unbiased error-phase weighted estimators. 



where 



Pi = 



100 
W 



J2 W 0) 



2 



and X(h) < X(k+i) Vfc < n. 



(Cll) 



APPENDIX D: 
ESTIMATORS 



SCATTER PLOTS OF RELATED 



Figures |D1[)D2| compare statistical parameters which pro- 
vide similar information, such as robust versus non-robust 
or normalized versus non-normalized estimators, as a func- 
tion of the weighting scheme. Error-phase weighted (and 
noise-unbiased, when applicable) estimators correspond to 
the most strongly peaked distributions around the correct 
values. 



APPENDIX E: UNWEIGHTED ESTIMATORS 
AS A FUNCTION OF VARIABILITY TYPE 

Figure |E1| presents a selection of unweighted estimators for 
stars of different variability types. This Figure is intended to 
be compared to Fig. f?\ which illustrates the same informa- 
tion for phase-error weighted and noise-unbiased estimators. 

This paper has been typeset from a TpX/ DTpX file prepared 
by the author. 
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Figure Dl. Deviations from the true values of (a) mean and median, (b) variance and interquartile range, (c) skewness and (d) kurtosis 
moments standardized by the estimated and true variances. The triangles and histograms in green denote error-weighted estimators, the 
crosses and histograms in red indicate error-phase weighted estimators, and the circles and histograms in blue represent noise-unbiased 
error-phase weighted estimators. 
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Figure El. A selection of unweighted estimators employed for classification is illustrated as a function of variability type. Class labels 
are described in Table [l] and denoted by symbols as shown in the legend of each panel. 
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